Week 4: Unary and Binary Operations

PPOL 6805 / DSAN 6750: GIS for Spatial Data Science
Fall 2025

Class Sessions

Author

Affiliation

Jeff Jacobs

jj1088@georgetown.edu

Published

Tuesday, September 17, 2024

Open slides in new tab →

Coordinate Reference Systems

First things first… what are coordinates?
We need a way to unambiguously specify where things are on the earth

Latitude and Longitude are Angles!

Angular Distance vs. Travel Distance

The Earth’s “width” is slightly greater than its “length” 😰

Projections

Smooshing 3D into 2D

Avoid Getting Lost in the Sauce

How To Avoid Getting Lost in the Sauce

**Tissot Circles**: Imagine infinitely small ellipses placed at regular intervals on the curved surface of the earth. Imagine these ellipses being projected along with the earth’s surface. When scaled up, changes in the ellipses show the location and quality of distortions on the projected map.

Choropleths for Good and Evil

Having Domain Knowledge of a Region

Crime in Mongolia

Population of Mongolia

Exhibit A

Exhibit B

Geospatial Operations 1: Unary Operations

Getting the Geometries

Using rnaturalearth with mapview

Code

set.seed(6805)
library(tidyverse)

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ lubridate 1.9.4     ✔ tibble    3.3.0
✔ purrr     1.0.4     ✔ tidyr     1.3.1
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Code

library(sf)

Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE

Code

library(rnaturalearth)
library(mapview)
france_sf <- ne_countries(country = "France", scale = 50)
(france_map <- mapview(france_sf, label = "geounit", legend = FALSE))

Centroid of France

Code

france_cent_sf <- sf::st_centroid(france_sf)

Warning: st_centroid assumes attributes are constant over geometries

Code

france_map + mapview(france_cent_sf, label = "Centroid", legend = FALSE)

One We Already Saw: Union

Computing the union of all geometries in the sf via sf::st_union()

Code

library(leaflet.extras2)

Loading required package: leaflet

Code

africa_sf <- ne_countries(continent = "Africa", scale = 50)
africa_union_sf <- sf::st_union(africa_sf)
africa_map <- mapview(africa_sf, label="geounit", legend=FALSE)
africa_union_map <- mapview(africa_union_sf, label="st_union(africa)", legend=FALSE)
africa_map | africa_union_map

Helpful for Rasterizing: BBox

Code

africa_bbox_sf <- sf::st_bbox(africa_sf)
africa_bbox_map <- mapview(africa_bbox_sf, label="st_bbox(africa)", legend=FALSE)
africa_map | africa_bbox_map

Convex Hulls by Country

Code

africa_countries_cvx <- sf::st_convex_hull(africa_sf)
africa_countries_cvx_map <- mapview(africa_countries_cvx, label="geounit", legend=FALSE)
africa_map | africa_countries_cvx_map

Convex Hull of Continent

Use st_union() first:

Code

africa_cvx <- africa_sf |> st_union() |> st_convex_hull()
africa_cvx_map <- mapview(africa_cvx, label="geounit", legend=FALSE)
africa_map | africa_cvx_map

One We Already Saw: Centroids

Computing the centroid of all geometries in the sf via sf::st_centroid()

Code

africa_cents_sf <- sf::st_centroid(africa_sf)

Warning: st_centroid assumes attributes are constant over geometries

Code

africa_cents_map <- mapview(africa_cents_sf, label="geounit", legend=FALSE)
africa_map | africa_cents_map

Geospatial Operations 2: Binary Operations

Spatial Joins

Code

nc <- system.file("shape/nc.shp", package="sf") |>
  read_sf() |>
  st_transform('EPSG:2264')
gr <- st_sf(
         label = apply(expand.grid(1:10, LETTERS[10:1])[,2:1], 1, paste0, collapse = ""),
         geom = st_make_grid(nc))
gr$col <- sf.colors(10, categorical = TRUE, alpha = .3)
# cut, to verify that NA's work out:
gr <- gr[-(1:30),]
suppressWarnings(nc_j <- st_join(nc, gr, largest = TRUE))
par(mfrow = c(2,1), mar = rep(0,4))
plot(st_geometry(nc_j), border = 'grey')
plot(st_geometry(gr), add = TRUE, col = gr$col)
text(st_coordinates(st_centroid(st_geometry(gr))), labels = gr$label, cex = .85)
# the joined dataset:
plot(st_geometry(nc_j), border = 'grey', col = nc_j$col)
text(st_coordinates(st_centroid(st_geometry(nc_j))), labels = nc_j$label, cex = .7)
plot(st_geometry(gr), border = '#88ff88aa', add = TRUE)

Spatial Sampling

Code

# Sample random points
africa_points_list <- sf::st_sample(africa_union_sf, 10)
africa_points_sf <- sf::st_sf(africa_points_list)
africa_points_map <- mapview(africa_points_sf, label="Random Point", col.regions=cb_palette[1], legend=FALSE)

Warning in cbind(`Feature ID` = fid, mat): number of rows of result is not a
multiple of vector length (arg 1)

Code

africa_map + africa_points_map

The “Default” Predicate: `st_intersects`

Code

countries_w_points <- africa_sf[africa_points_sf,]
mapview(countries_w_points, label="geounit", legend=FALSE) + africa_points_map

Counting with `lengths()`

Code

country_inter <- sf::st_intersects(africa_sf, africa_points_sf)
# Computes point counts for each polygon
(num_intersections <- lengths(country_inter))

 [1] 0 0 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 0 0 0
[39] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 2 2

Code

africa_sf <- africa_sf |> mutate(
  num_points = num_intersections
) |> arrange(geounit)
africa_sf |> select(geounit, num_points) |> head()

geounit	num_points	geometry
Algeria	2	MULTIPOLYGON (((8.576563 36…
Angola	2	MULTIPOLYGON (((13.07275 -4…
Benin	0	MULTIPOLYGON (((1.622656 6….
Botswana	0	MULTIPOLYGON (((25.25879 -1…
Burkina Faso	0	MULTIPOLYGON (((0.9004883 1…
Burundi	0	MULTIPOLYGON (((30.55361 -2…

Plotting with `mapview`

Code

mapview(africa_sf, zcol="num_points")

Plotting with `ggplot2`

Since we’re starting to get into data attributes rather than geometric features, switching to ggplot2 is recommended!

Code

africa_sf |> ggplot(aes(fill=num_points)) +
  geom_sf() +
  theme_classic()

One More Unary Operation: `st_buffer()`

Think about how you might answer questions like:
- “How far is your house (POINT) from Manhattan (POLYGON)?”
- “Are there any chemical plants within a mile of this building (POLYGON) / stretch of road (LINESTRING)?”
Lazy mode (my favorite mode): Compute distances from the centroid
GIS master mode: Construct the buffer!

On `POLYGON`s

Key line: manhattan_buffer_sf <- manhattan_union_sf |> st_buffer(dist = 1609.34) (1 mile \(\approx\) 1609.34 meters)

Code

library(tidyverse)
library(sf)
library(tidycensus)
library(tigris)
library(mapview)
options(tigris_use_cache = TRUE)
manhattan_sf <- get_acs(
  geography = "tract",
  variables = "B19013_001",
  state = "NY",
  county = "New York",
  year = 2020,
  geometry = TRUE,
  cb = FALSE
)
# Erase the island tracts real quick
island_tracts <- c(
  "Census Tract 1, New York County, New York",
  "Census Tract 2.02, New York County, New York"
)
manhattan_sf <- manhattan_sf |> filter(
  !(NAME %in% island_tracts)
)
# Union of all census tracts within the county
manhattan_union_sf <- st_union(manhattan_sf)
manhattan_union_map <- mapview(manhattan_union_sf, label="New York County")
# Construct buffer (1 mile ~= 1609.34 meters)
manhattan_buffer_sf <- manhattan_union_sf |> st_buffer(dist = 1609.34)
manhattan_buffer_map <- mapview(manhattan_buffer_sf, label="Buffer (1 Mile)", col.regions = cbPalette[1], legend = TRUE)
manhattan_buffer_map + manhattan_union_map

On `LINESTRING`s

Code

library(tidyverse)
library(sf)
## st_buffer, style options (taken from rgeos gBuffer)
l1 = st_as_sfc("LINESTRING(0 0,1 5,4 5,5 2,8 2,9 4,4 6.5)")
op = par(mfrow=c(2,3))
plot(st_buffer(l1, dist = 1, endCapStyle="ROUND"), reset = FALSE, main = "endCapStyle: ROUND")
plot(l1,col='blue',add=TRUE)
plot(st_buffer(l1, dist = 1, endCapStyle="FLAT"), reset = FALSE, main = "endCapStyle: FLAT")
plot(l1,col='blue',add=TRUE)
plot(st_buffer(l1, dist = 1, endCapStyle="SQUARE"), reset = FALSE, main = "endCapStyle: SQUARE")
plot(l1,col='blue',add=TRUE)
plot(st_buffer(l1, dist = 1, nQuadSegs=1), reset = FALSE, main = "nQuadSegs: 1")
plot(l1,col='blue',add=TRUE)
plot(st_buffer(l1, dist = 1, nQuadSegs=2), reset = FALSE, main = "nQuadSegs: 2")
plot(l1,col='blue',add=TRUE)
plot(st_buffer(l1, dist = 1, nQuadSegs= 5), reset = FALSE, main = "nQuadSegs: 5")
plot(l1,col='blue',add=TRUE)

What Makes Binary Operations “Fancier”?

Unary

st_centroid()
- POLYGON \(\mapsto\) POINT
- MULTIPOLYGON \(\mapsto\) POINT
st_convex_hull()
- POLYGON \(\mapsto\) POLYGON
- MULTIPOINT \(\mapsto\) POLYGON

Binary

st_intersection()
- (POINT, POINT) \(\mapsto\) POINT | POINT EMPTY
- (POLYGON, POLYGON) \(\mapsto\) POLYGON | LINESTRING | POINT | POLYGON EMPTY

st_is_empty() and st_dimension() become your new best friends 😉
st_is_empty(): Distinguishes between, e.g., POINT EMPTY and POINT(0 0)
st_dimension(): NA for empty versions, otherwise
- 2 for surfaces (POLYGON, MULTIPOLYGON)
- 1 for lines (LINESTRING, MULTILINESTRING)
- 0 for points (POINT, MULTIPOINT)

The Bad Kind of Overthinking: Will My Life Just Get Harder and Harder?

Unary Operations

Binary Operations

Ternary Operations

Quaternary Operations

Good News and Bad News

The good news: No!
The bad news: You’ll have to read the 465-page Volume I and then the 451-page Volume II and then to page 15 of Volume III of Cohn (1965) to know why:

(i spent 4 yrs of undergrad studying abstract algebra and now it all sits gathering dust somewhere deep within my brain plz just let me have this moment i’ll never mention it again i promise)

The Good Kind of Overthinking…

For fancier geospatial operations, we’ll need to start overthinking, about the possible relationships between two (or more) geometries! \(\leadsto\) Relational Predicates:

Figure 4.2 in Lovelace, Nowosad, and Muenchow (2024)

Ad Hoc / Ambiguous \(\rightarrow\) Precise: Enter DE-9IM!

DE-9IM Strings

Each cell here visualizes one component of the DE-9IM string 1020F1102, which describes the relationship between the following two geometries:

Boxey McBoxface: POLYGON(0 0, 1 0, 1 1, 0 1, 0 0)
Liney McLineface: LINESTRING(0.5 -0.5, 0.5 0.5)

Code

library(sf)
polygon <- po <- st_polygon(list(rbind(c(0,0), c(1,0), c(1,1), c(0,1), c(0,0))))
p0 <- st_polygon(list(rbind(c(-1,-1), c(2,-1), c(2,2), c(-1,2), c(-1,-1))))
line <- li <- st_linestring(rbind(c(.5, -.5), c(.5, 0.5)))
s <- st_sfc(po, li)

par(mfrow = c(3,3))
par(mar = c(1,1,1,1))

# "1020F1102"
# 1: 1
plot(s, col = c(NA, 'darkgreen'), border = 'blue', main = expression(paste("I(pol)",intersect(),"I(line) = 1")))
lines(rbind(c(.5,0), c(.5,.495)), col = 'red', lwd = 2)
points(0.5, 0.5, pch = 1)

# 2: 0
plot(s, col = c(NA, 'darkgreen'), border = 'blue', main = expression(paste("I(pol)",intersect(),"B(line) = 0")))
points(0.5, 0.5, col = 'red', pch = 16)

# 3: 2
plot(s, col = c(NA, 'darkgreen'), border = 'blue', main = expression(paste("I(pol)",intersect(),"E(line) = 2")))
plot(po, col = '#ff8888', add = TRUE)
plot(s, col = c(NA, 'darkgreen'), border = 'blue', add = TRUE)

# 4: 0
plot(s, col = c(NA, 'darkgreen'), border = 'blue', main = expression(paste("B(pol)",intersect(),"I(line) = 0")))
points(.5, 0, col = 'red', pch = 16)

# 5: F
plot(s, col = c(NA, 'darkgreen'), border = 'blue', main = expression(paste("B(pol)",intersect(),"B(line) = F")))

# 6: 1
plot(s, col = c(NA, 'darkgreen'), border = 'blue', main = expression(paste("B(pol)",intersect(),"E(line) = 1")))
plot(po, border = 'red', col = NA, add = TRUE, lwd = 2)

# 7: 1
plot(s, col = c(NA, 'darkgreen'), border = 'blue', main = expression(paste("E(pol)",intersect(),"I(line) = 1")))
lines(rbind(c(.5, -.5), c(.5, 0)), col = 'red', lwd = 2)

# 8: 0
plot(s, col = c(NA, 'darkgreen'), border = 'blue', main = expression(paste("E(pol)",intersect(),"B(line) = 0")))
points(.5, -.5, col = 'red', pch = 16)

# 9: 2
plot(s, col = c(NA, 'darkgreen'), border = 'blue', main = expression(paste("E(pol)",intersect(),"E(line) = 2")))
plot(p0 / po, col = '#ff8888', add = TRUE)
plot(s, col = c(NA, 'darkgreen'), border = 'blue', add = TRUE)

The predicate equals corresponds to the DE-9IM string "T*F**FFF*". If any two geometries obey this relationship, they are (topologically) equal!

Slowing Down: 9IM (no DE yet!)

Dimensionally Extended (DE) 9IM

Crunching it Down into a Tiny Box

DE-9IM	Interior	Boundary	Exterior
Interior	2	1	2
Boundary	1	0	1
Exterior	2	1	2

And Then into a Tiny String

212101212

And Then into an Infinitesimally-Small Point

DE-9IM Masks

Now terms can be given unambiguous, precise meaning!

`st_overlaps()`	Interior	Boundary	Exterior
Interior	`T`	`*`	`T`
Boundary	`*`	`*`	`*`
Exterior	`T`	`*`	`*`

Special Values (besides 0, 1, 2):
- T: “True” (non-empty, st_dimension() >= 0)
- F: “False” (empty, st_dimension() == NA)
- *: “Wildcard” (Don’t care what the value is)
st_overlaps(): T*T***T**, st_equals(): T*F**FFF*

DE-9IM vs. Everyday Language

DE-9IM can (in theory) represent \(6^9 = 10~077~696\) possible geometric relationships!
The English language has like 10, and they’re ambiguous ☠️ (Compromise employed by GIS systems: allow multiple masks for same English word:)

English	Mask	`212101212`	Result
“Disjoint”	`FFFF***`	`FALSE`	`x` not disjoint from `y`
“Touches”	`FT*******`	`FALSE`	`x` doesn’t touch `y`
“Touches”	`F*T**`	`FALSE`	`x` doesn’t touch `y`
“Crosses”	`TTT`	`TRUE`	`x` crosses `y`
“Within”	`TFF****`	`FALSE`	`x` is not within `y`
“Overlaps”	`TTT`	`TRUE`	`x` overlaps `y`

`st_relate()`: The Ultimate Predicate

Code

library(tidyverse)
library(rnaturalearth)
us <- ne_states(country="United States of America")
dc <- us |> filter(iso_3166_2 == "US-DC")
us <- us |>
  mutate(
    de9im = st_relate(us, dc),
    touch = st_touches(us, dc, sparse = F)
  ) |>
  select(iso_3166_2, name, de9im, touch) |>
  arrange(name)

although coordinates are longitude/latitude, st_relate assumes that they are
planar

Code

us

iso_3166_2	name	de9im	touch	geometry
US-AL	Alabama	FF2FF1212	FALSE	MULTIPOLYGON (((-87.41958 3…
US-AK	Alaska	FF2FF1212	FALSE	MULTIPOLYGON (((-141.0056 6…
US-AZ	Arizona	FF2FF1212	FALSE	MULTIPOLYGON (((-111.0063 3…
US-AR	Arkansas	FF2FF1212	FALSE	MULTIPOLYGON (((-90.30422 3…
US-CA	California	FF2FF1212	FALSE	MULTIPOLYGON (((-114.7243 3…
US-CO	Colorado	FF2FF1212	FALSE	MULTIPOLYGON (((-109.0463 4…
US-CT	Connecticut	FF2FF1212	FALSE	MULTIPOLYGON (((-73.6417 41…
US-DE	Delaware	FF2FF1212	FALSE	MULTIPOLYGON (((-75.05809 3…
US-DC	District of Columbia	2FFF1FFF2	FALSE	MULTIPOLYGON (((-77.02293 3…
US-FL	Florida	FF2FF1212	FALSE	MULTIPOLYGON (((-87.44734 3…
US-GA	Georgia	FF2FF1212	FALSE	MULTIPOLYGON (((-80.89029 3…
US-HI	Hawaii	FF2FF1212	FALSE	MULTIPOLYGON (((-154.8996 1…
US-ID	Idaho	FF2FF1212	FALSE	MULTIPOLYGON (((-117.0382 4…
US-IL	Illinois	FF2FF1212	FALSE	MULTIPOLYGON (((-89.1237 36…
US-IN	Indiana	FF2FF1212	FALSE	MULTIPOLYGON (((-84.80608 4…
US-IA	Iowa	FF2FF1212	FALSE	MULTIPOLYGON (((-96.48266 4…
US-KS	Kansas	FF2FF1212	FALSE	MULTIPOLYGON (((-102.0396 3…
US-KY	Kentucky	FF2FF1212	FALSE	MULTIPOLYGON (((-89.42446 3…
US-LA	Louisiana	FF2FF1212	FALSE	MULTIPOLYGON (((-89.52599 3…
US-ME	Maine	FF2FF1212	FALSE	MULTIPOLYGON (((-71.08495 4…
US-MD	Maryland	FF2F11212	TRUE	MULTIPOLYGON (((-75.64786 3…
US-MA	Massachusetts	FF2FF1212	FALSE	MULTIPOLYGON (((-71.19396 4…
US-MI	Michigan	FF2FF1212	FALSE	MULTIPOLYGON (((-84.4913 46…
US-MN	Minnesota	FF2FF1212	FALSE	MULTIPOLYGON (((-97.22609 4…
US-MS	Mississippi	FF2FF1212	FALSE	MULTIPOLYGON (((-88.40221 3…
US-MO	Missouri	FF2FF1212	FALSE	MULTIPOLYGON (((-95.31725 4…
US-MT	Montana	FF2FF1212	FALSE	MULTIPOLYGON (((-116.0482 4…
US-NE	Nebraska	FF2FF1212	FALSE	MULTIPOLYGON (((-104.0537 4…
US-NV	Nevada	FF2FF1212	FALSE	MULTIPOLYGON (((-114.0425 4…
US-NH	New Hampshire	FF2FF1212	FALSE	MULTIPOLYGON (((-71.50585 4…
US-NJ	New Jersey	FF2FF1212	FALSE	MULTIPOLYGON (((-75.54133 3…
US-NM	New Mexico	FF2FF1212	FALSE	MULTIPOLYGON (((-108.1375 3…
US-NY	New York	FF2FF1212	FALSE	MULTIPOLYGON (((-79.06523 4…
US-NC	North Carolina	FF2FF1212	FALSE	MULTIPOLYGON (((-76.03173 3…
US-ND	North Dakota	FF2FF1212	FALSE	MULTIPOLYGON (((-104.0476 4…
US-OH	Ohio	FF2FF1212	FALSE	MULTIPOLYGON (((-80.52023 4…
US-OK	Oklahoma	FF2FF1212	FALSE	MULTIPOLYGON (((-103.0002 3…
US-OR	Oregon	FF2FF1212	FALSE	MULTIPOLYGON (((-124.4924 4…
US-PA	Pennsylvania	FF2FF1212	FALSE	MULTIPOLYGON (((-79.76301 4…
US-RI	Rhode Island	FF2FF1212	FALSE	MULTIPOLYGON (((-71.23686 4…
US-SC	South Carolina	FF2FF1212	FALSE	MULTIPOLYGON (((-78.57316 3…
US-SD	South Dakota	FF2FF1212	FALSE	MULTIPOLYGON (((-104.0567 4…
US-TN	Tennessee	FF2FF1212	FALSE	MULTIPOLYGON (((-90.30422 3…
US-TX	Texas	FF2FF1212	FALSE	MULTIPOLYGON (((-103.3115 2…
US-UT	Utah	FF2FF1212	FALSE	MULTIPOLYGON (((-111.0502 4…
US-VT	Vermont	FF2FF1212	FALSE	MULTIPOLYGON (((-73.35134 4…
US-VA	Virginia	FF2F11212	TRUE	MULTIPOLYGON (((-76.01325 3…
US-WA	Washington	FF2FF1212	FALSE	MULTIPOLYGON (((-122.753 48…
US-WV	West Virginia	FF2FF1212	FALSE	MULTIPOLYGON (((-82.58945 3…
US-WI	Wisconsin	FF2FF1212	FALSE	MULTIPOLYGON (((-87.80425 4…
US-WY	Wyoming	FF2FF1212	FALSE	MULTIPOLYGON (((-109.0463 4…

(If You Don’t Want to Scroll)

Code

us |> filter(touch)

Warning: Using one column matrices in `filter()` was deprecated in dplyr 1.1.0.
ℹ Please use one dimensional logical vectors instead.

iso_3166_2	name	de9im	touch	geometry
US-MD	Maryland	FF2F11212	TRUE	MULTIPOLYGON (((-75.64786 3…
US-VA	Virginia	FF2F11212	TRUE	MULTIPOLYGON (((-76.01325 3…

Doing Things with DE-9IM (Back to Binary Operations)

Almost a Spatial Join

Code

library(mapview)
library(leaflet.extras2)
africa_sf <- ne_countries(continent = "Africa", scale = 50) |> select(iso_a3, geounit)
africa_map <- mapview(africa_sf, label="geounit", legend=FALSE)
N <- 10
africa_union_sf <- sf::st_union(africa_sf)
sampled_points_sf <- sf::st_sample(africa_union_sf, N) |> sf::st_sf() |> mutate(temp = runif(N, 0, 100))
sampled_points_map <- mapview(sampled_points_sf, label="Random Point", col.regions=cbPalette[1], legend=FALSE)
countries_points_sf <- africa_sf[sampled_points_sf,]
filtered_map <- mapview(countries_points_sf, label="geounit", legend=FALSE) + sampled_points_map
(africa_map + sampled_points_map) | filtered_map

Spatial Filter \(\neq\) Spatial Join

The issue: Data attributes of POINTs are not merged into data attributes of POLYGONs

POINT Attributes

Code

st_geometry(sampled_points_sf) <- c("geom")
sampled_points_sf |> head()

geom	temp
POINT (-5.135802 22.30373)	99.14195
POINT (48.01641 10.55753)	41.21422
POINT (47.30147 5.038318)	95.30175
POINT (22.64591 7.799242)	68.37119
POINT (7.739425 5.95728)	56.40943
POINT (22.13408 20.87522)	86.12631

POLYGON Attributes

Code

countries_points_sf |> head(4)

	iso_a3	geounit	geometry
2	ZMB	Zambia	MULTIPOLYGON (((30.39609 -1…
58	SOM	Somalia	MULTIPOLYGON (((41.53271 -1…
59	-99	Somaliland	MULTIPOLYGON (((48.93857 11…
91	NGA	Nigeria	MULTIPOLYGON (((7.300781 4….

Our First Real Spatial Join: `st_join()`

Code

joined_sf <- countries_points_sf |> st_join(sampled_points_sf)
joined_sf |> head()

	iso_a3	geounit	temp	geometry
2	ZMB	Zambia	5.788125	MULTIPOLYGON (((30.39609 -1…
58	SOM	Somalia	95.301751	MULTIPOLYGON (((41.53271 -1…
59	-99	Somaliland	41.214222	MULTIPOLYGON (((48.93857 11…
91	NGA	Nigeria	56.409432	MULTIPOLYGON (((7.300781 4….
114	MLI	Mali	99.141948	MULTIPOLYGON (((-11.3894 12…
123	LBY	Libya	86.126311	MULTIPOLYGON (((9.51875 30….

But… We Were Still in Easy Mode

Every point could be matched one-to-one with a country. But what if… 😱

Code

g <- st_make_grid(st_bbox(st_as_sfc("LINESTRING(0 0,1 1)")), n = c(2,2))
par(mar = rep(0,4))
plot(g)
plot(g[1] * diag(c(3/4, 1)) + c(0.25, 0.125), add = TRUE, lty = 2)
text(c(.2, .8, .2, .8), c(.2, .2, .8, .8), c(1,2,4,8), col = 'red')

Spatially Intensive vs. Spatially Extensive

Extensive attributes: associated with a physical size (length, area, volume, counts of items). Ex: population count.
- Associated with an area \(\implies\) if that area is cut into smaller areas, the population count needs to be split too
- (At minimum, the sum of the population counts for the smaller areas needs to equal the total for the larger area)
Intensive attributes: Not proportional to support: if the area is split, values may vary but on average remain the same. Ex: population density
- If an area is split into smaller areas, population density is not split similarly!
- The sum of population densities for the smaller areas is a meaningless measure
- Instead, the mean will be more useful as ~similar to the density of the total

Handling the Extensive Case

Assume the extensive attribute \(Y\) is uniformly distributed over a space \(S_i\) (e.g., for population counts we assume everyone is evenly-spaced across the region)
We first compute \(Y_{ij}\), derived from \(Y_i\) for a sub-area of \(S_i\), \(A_{ij} = S_i \cap T_j\):

\[ \hat{Y}_{ij}(A_{ij}) = \frac{|A_{ij}|}{|S_i|}Y_i(S_i) \]

where \(|\cdot|\) denotes area.
Then we can compute \(Y_j(T_j)\) by summing all the elements over area \(T_j\):

\[ \hat{Y}_j(T_j) = \sum_{i=1}^{p}\frac{|A_{ij}|}{|S_i|}Y_i(S_i) \]

Handling the Intensive Case

Assume the variable \(Y\) has constant value over a space \(S_i\) (e.g., population density in assumed to be the same across all sub-areas)
Then the estimate for a sub-area is the same as the estimate for the total area:

\[ \hat{Y}_{ij} = Y_i(S_i) \]

So that we can obtain estimates of \(Y\) for new spatial units \(T_j\) via area-weighted average of the source values:

\[ \hat{Y}_j(T_j) = \sum_{i=1}^{p}\frac{|A_{ij}|}{|T_j|}Y_j(S_i) \]

References

Cohn, P. M. 1965. Universal Algebra. Springer Science & Business Media.

Krygier, John, and Denis Wood. 2016. Making Maps, Third Edition: A Visual Guide to Map Design for GIS. Guilford Publications.

Lovelace, Robin, Jakub Nowosad, and Jannes Muenchow. 2024. Geocomputation with R. CRC Press. https://r.geocompx.org/.

Monmonier, Mark. 2018. How to Lie with Maps. University of Chicago Press. https://www.dropbox.com/scl/fi/7rsqbxgge6llggnaf5tit/How-to-Lie-with-Maps-Third-Edition.pdf?rlkey=6rqxta7cjyq3oqdnskj4dtwj8&dl=1.

Pebesma, Edzer, and Roger Bivand. 2023. Spatial Data Science: With Applications in R. CRC Press. https://r-spatial.org/book/.

Coordinate Reference Systems

Latitude and Longitude are Angles!

Angular Distance vs. Travel Distance

Projections

Smooshing 3D into 2D

Avoid Getting Lost in the Sauce

How To Avoid Getting Lost in the Sauce

Choropleths for Good and Evil

Having Domain Knowledge of a Region

Crime in Mongolia

Population of Mongolia

Exhibit A

Exhibit B

Geospatial Operations 1: Unary Operations

Getting the Geometries

Centroid of France

One We Already Saw: Union

Helpful for Rasterizing: BBox

Convex Hulls by Country

Convex Hull of Continent

One We Already Saw: Centroids

Geospatial Operations 2: Binary Operations

Spatial Joins

Spatial Sampling

The “Default” Predicate: st_intersects

Counting with lengths()

Plotting with mapview

Plotting with ggplot2

One More Unary Operation: st_buffer()

On POLYGONs

On LINESTRINGs

What Makes Binary Operations “Fancier”?

The Bad Kind of Overthinking: Will My Life Just Get Harder and Harder?

Good News and Bad News

The Good Kind of Overthinking…

Ad Hoc / Ambiguous \(\rightarrow\) Precise: Enter DE-9IM!

DE-9IM Strings

Slowing Down: 9IM (no DE yet!)

Dimensionally Extended (DE) 9IM

Crunching it Down into a Tiny Box

And Then into a Tiny String

And Then into an Infinitesimally-Small Point

DE-9IM Masks

DE-9IM vs. Everyday Language

st_relate(): The Ultimate Predicate

(If You Don’t Want to Scroll)

Doing Things with DE-9IM (Back to Binary Operations)

Almost a Spatial Join

Spatial Filter \(\neq\) Spatial Join

Our First Real Spatial Join: st_join()

But… We Were Still in Easy Mode

Spatially Intensive vs. Spatially Extensive

Handling the Extensive Case

Handling the Intensive Case

References

The “Default” Predicate: `st_intersects`

Counting with `lengths()`

Plotting with `mapview`

Plotting with `ggplot2`

One More Unary Operation: `st_buffer()`

On `POLYGON`s

On `LINESTRING`s

`st_relate()`: The Ultimate Predicate

Our First Real Spatial Join: `st_join()`